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Implementation of Z-vector method in the relativistic coupled cluster framework to 
calculate first order energy derivatives: Application to SrF molecule 
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The molecular dipole moment and magnetic hyperfine structure constant demand an accurate 
wavefunction far from the nucleus and in near nuclear region, respectively. We, therefore, employ 
the so-called Z-vector method in the domain of relativistic coupled cluster theory to calculate the first 
order property of molecular systems in their open-shell ground state configuration. The implemented 
method is applied to calculate molecular dipole moment and parallel component of the magnetic 
hyperfine structure constant of SrF molecule. The results of our calculation are compared with the 
experimental and other available theoretically calculated values. We are successful in achieving good 
accordance with the experimental results. The result of our calculation of molecular dipole moment 
is in the accuracy of ~ 0.5 %, which is clearly an improvement over the previous calculation based 
on the expectation value method in the four component coupled cluster framework [V. S. Prasannaa 
et al, Phys. Rev. A 90, 052507 (2014)] and it is the best calculated value till date. Thus, it can 
be inferred that the Z-vector method can provide an accurate wavefunction in both near and far 
nuclear region, which is evident from our calculated results. 
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Theoretical physicists find it very challenging to calcu¬ 
late the spectroscopic properties of atoms and molecules. 
The precise description of the spectroscopic properties 
demands the wavefunction to be accurate both in the 
nuclear region and the region far from the nucleus. The 
calculation of an accurate wavefunction involving heavy 
atoms and molecules needs to include the relativistic and 
electron correlation effects simultaneously, as these two 
effects are non-additive in nature [2,0/ The best pos¬ 
sible way to include the effects of relativity in a single 
determinantal theory is to solve the Dirac-Hartree-Fock 
(DHF) Hamiltonian in its four component formalism. 
The DHF Hamiltonian converts the complicated many 
electron problem into a sum of many one-electron prob¬ 
lems by assuming an average electron-electron interac¬ 
tion. Therefore, the DHF Hamiltonian lacks the correla¬ 
tion of opposite spin electrons. The missing electron cor¬ 
relation can be included by adding orthogonal space to 
the DHF wavefunction. On the other hand, normal cou¬ 
pled cluster (NCC) @-0 method is known to be the most 
elegant many-body theory to effectuate the dynamic part 
of the electron correlation. 

The calculations of one electron response properties in 
the NCC framework, can either be done by taking ex¬ 
pectation value of the desired property operator or as a 
derivative of energy. These two approaches are not same 
as the NCC is by nature non-variational. In fact, the 
first order derivative of energy is the corresponding ex¬ 
pectation value plus some additional terms, which makes 
the derivative approach closer to the full configuration 
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interaction (FCI) property value. It is worth to mention 
that the expectation value approach in the NCC leads 
to a nonterminating series and any truncation scheme 
introduces an additional error [6j. 

In general, the energy is a function of both the deter¬ 
minantal coefficients (Co) in the expansion of the many 
electron correlated wavefunction and the molecular or¬ 
bital coefficients (Cm) for a fixed nuclear geometry. The 
first order energy derivative in NCC can be written as 

5E[C D (X),C M (X)\ SE 6C d 5E SC m 
SX " SC D SX + SC M SX ' 

Thus, for the calculation of energy derivative in NCC 
framework, it is, therefore, necessary to calculate the 
derivative of energy with respect the determinantal co¬ 
efficients as well as the molecular orbital coefficients. It 
further requires the derivative of the determinantal co¬ 
efficients and molecular orbital coefficients with respect 
to the external field of perturbation. However, Bartlett 
and co-workers 0 have shown that these derivative terms 
can be transformed into a single linear equation by us- 
ingZ-vector method. The advantage of Z-vector method 
[Si 0 is that for the calculation of several properties, one 
needs to solve a single linear equation instead of solving 
equations for each external perturbation field of interest. 
The detailed diagrammatic of Z-vector method in NCC 
framework is given in Ref. [§|. 

In this rapid communication, we have shown that the 
Z-vector method in the NCC framework within its four 
component description can generate an accurate wave- 
function in the near nuclear region as well as in the region 
far from the nucleus. To justify our argument, we have 
compared magnetic hyperfine structure constant (HFS) 
and molecular dipole moment of SrF with the experi- 
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mental values as the calculation of these properties need 
an accurate wavefunction in the near nuclear region and 
the region far from the nucleus, respectively. The rea¬ 
son for choosing SrF molecule is as follows: The knowl¬ 
edge of long range dipole-dipole interaction is very im- 
rortant to produce ultracold molecules in optical lattice 
K|. SrF molecule can be cooled by laser spectroscopy 
11] and thus it can be used for high precession spec¬ 
troscopy EH, Q- Currently an experimental search for 
parity violation using SrF is in progress 0- Therefore, 
detailed knowledge of the spectroscopic properties like 
dipole moment and magnetic HFS is very important to 
interpret the experimental findings. 

The Many body time-independent Dirac-Coulomb 
Hamiltonian is given by 

h = y, [ - ecu • v, + os - i)c 2 +vr c (n)+J2 — 

i A j>i 7 ' i? 

where c is the speed of light, a and /3 are the usual Dirac 
matrices, I is the 4x4 identity matrix. uc (ri) is the 
nuclear potential term of nucleus A. 

The dynamic part of the electron correlation is in¬ 
cluded using the coupled cluster method. The wave func¬ 
tion in the coupled cluster method is defined as 

l^cc) = e T |<f> 0 ) (2) 

where, | < f > o) is the ground state single determinant wave- 
function and T is cluster operator which is given by 

N 

T = T 1 + T 2 + ...+T N = J2 T n (3) 

n 

with 

Tm = (ml) 2 5Z i ij.'.' a a a b- — a 3 a i W 

' ' ij..ab.. 

i,j(a,b) are the hole (particle) indices and i“ b ;- are the clus¬ 
ter amplitudes corresponding to the cluster operator T m . 
The equations for n-body cluster amplitudes and corre¬ 
lation energy are given by 

($“;;|(H J ve T ) c |$ 0 )=0 (5) 

(<f> 0 \(H N e T )c\$o) = E corr ( 6 ) 

where Hn is the normal ordered Hamiltonian and the 

subscript C indicates only the connected terms in the 
contraction between Hjy and T. The connectedness en¬ 
sures size-extensivity. Once the cluster amplitude equa¬ 
tion © is solved, the correlation energy can be obtained 
from equation [Gj The coupled cluster energy is a func¬ 
tion of both the determinantal coefficients ( Cd) and the 
molecular orbital coefficients {Cm)- Therefore, the cal¬ 
culation of coupled cluster energy derivative need both 
the derivative of Cps and Cms with respect to external 


field of perturbation. However, the equations involving 
derivative of Cps and Cms are linear equations. Thus, 
one needs to solve the linear equations for each external 
field perturbation of interest. This can be avoided with 
the introduction of an antisymmetrized de-excitation op¬ 
erator, A, where the equations for the amplitudes of A 
is perturbation independent. Therefore, the solution of 
one linear equation is required instead of solving for each 
external perturbation. The second quantized form of the 
perturbation independent operator, A is given by 

N 

A = Ai + A 2 + ... + A at = A n (7) 

n 

where 

) = ( m !) 2 ^ab.. a i a j-"- a b a a ( 8 ) 

where i,j(a,b) are the hole(particle) indices and X^ h are 
the cluster amplitudes corresponding to the cluster oper¬ 
ator A m . The detailed description of A operator and A 
amplitude equation is given in reference @. The working 
A amplitude equation is given by 

($ 0 | [A(iJjve T )c]c|^“'.'> + ($o\{H N e T )c\<t>rnt) 

<^int | A|3>“") + ($o\(H N e T )c\$?::) = o (9) 

where is the determinant corresponding to the in¬ 
termediate excitation between $o and <3?“\ In the cou¬ 
pled cluster single and double (CCSD) model, A becomes, 
A = Ai + A 2 . The explicit equations for the amplitudes 
of Ai and A 2 operators are 

($o| [A{H N e T ) c ]cm) + ($o\{H N e T ) c m) = 0,(10) 

($o| [A{H N e T ) c ]c\<f>%) + (<S>o\{H N e T ) c m) 

($“|A|$“ b ) + (<S> 0 \{H N e T ) c \<t>%) = 0. (11) 

It is interesting to note that the term 
(4>o|(i?Are T )c|$“)($-|A|$“ b ) of equation ECU yields 
one disconnected diagram, which is given in Fig. [I] The 
said diagram is not of the type of closed with discon¬ 
nected part. This ensures that the energy derivative is 



FIG. 1: Disconnected yet linked diagram in A 2 equation. 
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linked and thus size extensive. The equation for energy 
derivative can be written as 


A E' = <$o|(Chve T ) c |$o) + <$o|[A(Oe T ) c ]c|$o) (12) 

where On is the derivative of normal ordered perturbed 
Hamiltonian with respect to external field of perturba¬ 
tion. 

The molecular dipole moment of a heavy diatomic 
molecule arises due to the fact that the nuclear charge 
of the two atoms are not same and thus the electron den¬ 
sity is not evenly distributed around the two nucleus. In 
the Born-Oppenheimer approximation, we can separate 
out the nuclear and electronic contribution from the total 
contribution. The dipole moment operator is given by 

P = - ^2 U + ^2 z a?a (13) 

i A 


where i stands for the electron and Za is the atomic num¬ 
ber of nucleus A. The first term of the above equation 
is the electronic contribution and the second term is the 
nuclear contribution. 

The magnetic HFS arises due to the interaction of nu¬ 
clear magnetic dipole moment with the magnetic moment 
of electrons. Thus, the magnetic HFS is the sum of one 
body interactions from the point of view of electronic 
structure theory @|. The magnetic vector potential (A) 
at a distance r is given by 


PkXr 

• 2 * ^3 


(14) 


where pk is the magnetic moment of the nucleus K. The 
perturbed HFS Hamiltonian of an atom due to A in the 
Dirac theory is given by Hh yv = Yl'i Pi'Aii where n is the 
total no of electrons and on denotes the Dirac a matrices 
for the \ th electron. The z projection (along the molec¬ 
ular axis) of the expectation value of the corresponding 
perturbed HFS Hamiltonian gives the parallel magnetic 
HFS constant, , which is given by 




pk 

in 




an x n 


I’M. 


(15) 


where I is the nuclear spin quantum number and H repre¬ 
sents the z component (along molecular axis) of the total 
angular momentum of the diatomic molecule. 

The N electron ground and excited determinants are 
constructed with the one electron spinors, those are the 
solutions of Dirac-Hartree-Fock equation. The DIRAC10 
0 program package is used to solve the Dirac-Fock 
equation and to obtain the matrix elements required for 
property calculations. Gaussian charge distribution is 
considered as the nuclear model where the nuclear pa¬ 
rameters [lbj are taken as default values in DIRAC10. 
Large and small component basis functions are gener¬ 
ated by applying restricted kinetic balance (RKB) |17| in 
which basis functions are represented in scalar basis and 


unphysical solutions are removed by diagonalizing the 
free particle Hamiltonian. This generates the electronic 
and positronic solution in 1:1 manner. We have done five 
different sets of calculation using five different basis sets 
for Sr and F. These are cc-pVDZ, aug-cc-pVTZ, aug-cc- 
pCVTZ, d-aug-cc-pCVTZ and aug-cc-pCVQZ for F atom 
fl8| and dyall.v2z, dyall.v3z, dyall.cv3z, d-aug-dyall.cv3z 
and dyall.cv4z for Sr atom [131. Both large and small 
component basis are taken in uncontracted form. None 
of the electrons are frozen in our correlation calculation 
and the virtual orbitals whose energies are greater than 
a certain threshold are not considered as the high energy 
virtual orbitals contribute less in correlation calculation. 

We have taken the following strategies to code the Z- 
vector method. First, the one electron and two electron 
matrix elements are obtained from DIRAC10 package 
0. Then we have solved the NCC part i.e., the Tj 
and X 2 amplitude equations. This is followed by the con¬ 
struction of different types of H (H = ( He T )c )■ After 
that H vertices are contracted with one Ai or A 2 vertex 
to construct the Ai and A 2 amplitude equations. At the 
end Tj, Tj, Ai and A 2 amplitudes are contracted with 
property integrals to get corresponding property value. 
To solve the T\ and T 2 amplitudes and to construct the 
H , we have used a recursive intermediate factorization 
of diagrams as described by Bartlett and coworkers 0] . 
This saves enormous computational cost. To debug this 
code, we benchmarked our correlation energy with the 
results obtained from DIRAC10 with same basis, same 
convergence criteria and using same direct inversion in 
the iterative subspace (DIIS). We have achieved 7 to 8 
decimal place agreement with DIRAC10 for correlation 
energy independent of the choice of molecules as well as 
of the basis sets. The discrepancy beyond this limit could 
be due to the use of cutoff in storing of the intermediate 
diagrams or the use of different convergence algorithm. 
The H and H matrix elements are stored by setting a 
cutoff of 10 -12 to save storage requirement as the con¬ 
tribution of the two body matrix elements beyond that 
limit is negligible. The tolerance used for the convergence 
of both T and A amplitudes is 10~ 9 . We have used the 
experimental bond length (2.075 A) of SrF 01 in all the 
calculations. 

In Table [I] we present the sets of basis used in our cal¬ 
culations and each combination is denoted by an English 
alphabet letter. The fourth and fifth column of Table U 
represent the cutoff used and the number of spinor gen¬ 
erated using that cutoff for correlation calculation, re¬ 
spectively. We also compiled the correlation energy of 
SrF obtained from CCSD and second-order many body 
perturbation theory (MBPT(2)), which uses a first-order 
perturbed wavefunction. 

In Table [Hi we present the molecular dipole moment in 
Debye of the ground state of SrF molecule in five different 
basis sets. The experimental value 0] is also presented 
in the same table for comparison. It is clear from the 
table that with increase in the number of basis function 
the dipole moment converges towards the experimental 
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TABLE I: Cutoff used and correlation energy of the ground state of SrF in different basis sets 


Basis Cutoff Spinor Correlation Energy (a.u.) 


Name 

Sr 

F 

(a.u.) 


MBPT(2) 

CCSD 

A 

dyall.v2z 

cc-pVDZ 


298 

-1.528622377 

-1.474175309 

B 

dyall.v3z 

aug-cc-pVTZ 

500 

366 

-1.292096889 

-1.256789931 

C 

dyall.cv3z 

aug-cc-pCVTZ 

500 

436 

-1.838672666 

-1.775710294 

D 

d-aug-dyall.cv3z 

d-aug-cc-pCVTZ 

100 

596 

-1.681326043 

-1.621603448 

E 

dyall.cv4z 

aug-cc-pCVQZ 

50 

520 

-1.454995437 

-1.402836367 


value. This is expected as more basis functions generate 
more correlation space and thereby improve the dipole 
moment. In particular, our calculated dipole moment 
in basis E is very close (~ 0.5 %) to the experimental 
value. The results obtained for the dipole moment of 


TABLE II: Molecular dipole moment (in Debye) of the 
ground state of SrF 


Basis 

Z-vector 

Experiment [22] 

A 

3.0158 


B 

3.3898 


C 

3.4023 

3.4676(10) 

D 

3.4376 


E 

3.4504 



TABLE III: Comparison of molecular dipole moment of the 
ground state of SrF in different methods 

Method 

Reference 

Dipole (D) 

Ionic model 

Torring et al [25] 

3.67 

SCF 

Langhoff et al. [24] 

2.579 

CPF 

Langhoff et al. [24] 

3.199 

CISD 

Langhoff et al. [24] 

2.523 

EPM 

Mestdagh et al. [25] 

3.6 

HF (finite difference) 

Kobus et al. [26] 

2.5759 

CCSD 

Prasannaa et al. [27] 

3.41 

Z-vector 

This work(E) 

3.4504 

Expt. 

Ernst et al. [22] 

3.4676(10) 


the ground state of SrF by other methods and experiment 
are compiled in Table Hill The dipole moment of SrF was 
first calculated by Torring et al, [23| by using an ionic 
model and they got a value of 3.67 D. Langhoff et al, (24| 
performed the first ab initio calculation of the dipole mo¬ 
ment of SrF by using Slater type of basis function. They 
reported the dipole moment using three different meth¬ 
ods, i.e., self consistent field (SCF), configuration inter¬ 
action in single and double approximation (CISD) and 
the coupled pair function (CPF) method. Among them, 
CISD method is not size extensive while CPF is, thus 


CPF approach gives better agreement with experiment. 
However, Langhoff et al, did not consider the relativistic 
motion of electrons. Mestdagh et al, [25| used electro¬ 
static polarization model and got 3.6 D as a molecular 
dipole moment of SrF. Kobus et al, [26[ obtained a dipole 
moment of 2.5759 D by using finite difference method in 
the Hartree-Fock (HF) level. The first relativistic calcu¬ 
lation of dipole moment of SrF in the CCSD model was 
calculated by Prasannaa et al, [23] by taking expectation 
value of the corresponding operator. The expectation 
value framework leads to a connected yet nonterminat¬ 
ing series. Prasannaa et al, took only the linear terms 
in the property calculations using CCSD wavefunction 
and got 3.41 D as a result. Our four component Z-vector 
calculation give s a result of 3.4504 D by using dyall.cv4z 
basis for Sr jlfij] and aug-cc-pCVQZ basis for F |19j ] (ba¬ 
sis E) and this result shows the best agreement with the 
experiment so far. 

In Table HVl we present the parallel component of the 
magnetic HFS constant of 8 'Sr and 19 F of the ground 
state of SrF molecule. We also present the experimental 
value [13] of those in the same table for comparison. Our 
calculated result using Z-vector method show good agree¬ 
ment with the experimental result. The highest and low¬ 
est deviation from experimental values for parallel mag¬ 
netic HFS constant of 87 Sr atom are for the A (~ 45 
MHz) and C (~ 24 MHz) basis respectively. For the par¬ 
allel magnetic HFS constant of 19 F, the maximum and 
minimum deviation occur for the basis D (~ 10 MHz) 
and A (~ 4 MHz) basis. 

The calculated magnetic HFS constant values are in 
good agreement with the sophisticated experiment but 
the extent of accuracy is not so in comparison to that of 
the calculated dipole moment values. This could possibly 
be due to the fact that as we proceed from basis A to E, 
we have added extra Gaussian type orbitals (GTOs) of 
higher angular momentum in addition to lower angular 
momentum. As the higher angular momentum GTO 
shifts the electron density towards the outer region, the 
addition of higher angular momentum GTO improves 
the outer region much better than the inner region of the 
molecular wavefunction. This is why as we go from basis 
A to E, our molecular dipole moment value matches 
more closely than the magnetic HFS values with the 
experimental results. 
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TABLE IV: Parallel (Ay) magnetic hyperfine structure constant of the ground state of SrF in MHz 


Basis 


5-i 

CO 

t- 

00 


19p 

Z-vector 

Experiment [28] 

Z-vector 

Experiment [28] 

A 

546.08 


121.93 


B 

558.96 


118.70 


C 

566.62 

591(3) 

119.64 

126(3) 

D 

561.25 


116.35 


E 

559.65 


117.74 



We have successfully implemented Z-vector method in 
the relativistic NCC domain using 4-component wave 
function to calculate the first order energy derivatives. 
We applied this method to calculate the molecular dipole 
moment and parallel magnetic HFS constant of SrF 
molecule. The results from our calculations are in good 
agreement with the experimental values. Therefore, we 


can conclude that Z-vector method in the relativistic 
framework can produce an accurate wavefunction in the 
near nuclear region as well as far from the nucleus. 
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